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yA/ms. We independently develop a simulation code following the previous dynamical Monte Carlo simulation of the diffusive shock 
acceleration under the isotropic scattering law during the scattering process, and the same results are obtained. 
Methods. Since the same results test the validity of the dynamical Monte Carlo method for simulating a collisionless shock, we extend 
the simulation toward including an anisotropic scattering law for further developing this dynamical Monte Carlo simulation. Under 
this extended anisotropic scattering law, a Gaussian distribution function is used to describe the variation of scattering angles in the 
particle's local frame. 

Results. As a result, we obtain a series of different shock structures and evolutions in terms of the standard deviation values of the 
given Gaussian scattering angular distributions. 

Conclusions. We find that the total energy spectral index increases as the standard deviation value of the scattering angular distribution 
increases, but the subshock's energy spectral index decreases as the standard deviation value of the scattering angular distribution 
increases. 

Key words, acceleration of particles- shock wave-solar energetic particle 
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It is well known that the diffusive accel eration model has been 
popular for more than five decades since iFermil (|1949|) first pro- 
posed that cosmic rays could be produced via diffusive pro- 
cesses. Until now, diffusive shock acceleration (DSA) (i.e. first- 
order Fermi acceleration) has been extensively applied to many 
physical systems, such as shocks in the solar system, in ou r 
Galaxy, and throughout the Universe (lBaring||l997tlBellll2004 . 

It is also well known that the nonlinear interactions in 
plasma usually include such things as the turbulence of scatter- 
ing wave field, cosmic ray (CR) injection, and "back reaction" 
by CR pressure. These complex behaviors have held back com- 
prehensive understanding of the DSA and nonlinear DSA the- 
ory. Therefore, to study the properties of the acceleration pro- 
\ cess and dynamical behavior of the CR's "back reaction" on 
■ the background flow, choosing numerical si mulation methods 
has been a pr i mary a nd essential pr o blem. (lOstrowskil [l99H 
Kang& Jones 1 99lt iMalkovl Il997b lEUison. Baring & Jones 
— - ' 1996T iKnerr. Jokipii & Ellisonl Il996h IBerezhko et al 



computational grid for sca ttering calculation is don e by pa rticle- 
in-cell (PIC) techniques. IKnerr. Jokipii & Ellisonl (Il996l) suc- 
cessfully developed the dynamical Monte Carlo simulations 
for the Earth's bow shock with important results in the maxi- 
mum energetic particles achieving greater than lMeV acceler- 
ated b y the shock. Before the dynamical M onte Carlo simula- 
tion, lEUison. Mobius & P aschmann (1990) had developed sta- 
tionary Monte Carlo simulations with the result that the cutoff 
energy accelerated by the shock only r eached 100 keV owing to 
the li mitation on the size of bow shock. Barin g. Ellison & Jones I 
(1 19951) also used the stationary Monte Carlo method for simulat- 
ing the oblique interplanetary shocks, and the calculated results 
are a good fit to the observed data. There are many works tha t 
use the Monte Carlo met hod: lEUison. Baring & Jonesl (fl996). 
Ellis on & Doublel (120021) . IVladimirov. Ellison & Bvkovl 72006. 
120081) . and others. 
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I Shalchi. Li & Zankll20l6l: IBerezhko & Volkll2000l) . The 
main simulation methods are introduce d here, and a more de- 
tailed review can be seen in lKangl (1200 ll) . 

Monte Carlo method. In Monte Carlo simulations, the scat- 
tering processes between individual particles with the collective 
background flow are simulated around a one-dimensional paral- 
lel shock. The particle scattering process is based on a prescribed 
scattering law, and collecting moments based on the background 
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Hybrid simulation. The particles' motion equations are 
solved explicitly based on the electromagnetic field of the back- 
ground plasma. Since the proton mass is about 2000 times that 
of the electron mass, the total plasma is assumed to be a coupling 
of two components: one component (e.g. electrons) is treated as 
a massless fluid and the other componentfe.g. protons, ions) is 
treated as individual particles dLeroy et al.lll982l) . This method 
also employs the PIC techniques based on computational grids. 
However, the limited computational resources imply that the ex- 
tensive calculation of the electromagnetic field uses unrealistic 
parameters and are unable to follow the shock for long enough 
(IGiacalone etaDl 19931: iGiacalone & Jokipiill2Q09h . 

Two-fluid model. The two-fluid model uses the diffusive- 
convection equation, coupled with the gas dynamic equa- 
tions, to simulate the CR's acceleration as a gas component 
and an accelerated particle component (iDrury & Falld 1 19861: 
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iDorfil 119901: Uones & Kan j 1 1992b . Since the CR energy den- 
sity is solved instead of the particle's distribution function 
in this model, the simulation results are based on some as- 
sumptions, such as the CR's adiabatic index, averaged diffu- 
sive coefficient, and injection rate. The averaged diffusive coef- 
ficient needs to be inferred from the diffusive-advection equa- 
tion. The acceleration efficiency dependens on the assump- 
tion of the injection rate. Under the reasonable a priori mod- 
els of these free parameters, a lot of simulations have tested 
the acceleration efficiency and the shock structure, and both 
agree with those derived from the diffusion-advection method. 
And these semi-analytic solutions t hat have been extensively 
used c a n be seen in recent w o rks: iMalkov. Diamond & V olk 
(l2QQQh iBednarz & Qstrowskil (l200l[). iBlasil (|2002l 12004 . 
lAmato & Blasil (120051 l2006h. iBlasi. Amato & Capriolil (120071) . 

l5J. 



ICaprioli et al.1 (l2009ah . lCaprioli. Amato & Blasil (12010a 

Kinetic simulation. Within this full numerical simulation, 
the diffusion-convection equation for the distribution function 
is solved with a momentum-dependent dif fusion coefficient and 
a suitable assumption of injection rate (jKang & Jonesl 119911: 
lBerezhkoetaDll994t iBerezhko & Volkll200Qh . Unlike the two- 
fluid model, the kinetic model should not assume the CR's 
adiabatic index, in addition to using the momentum-dependent 
diffusion coefficient instead of the averaged diffusion coeffi- 
cient. Berezhko and collaborators have studied numerous DSA 
models for supernova remnants (SNRs) with the kinetic model 
and conclude that about 20 % of the total SN energy trans- 
ferred to the CR's population in the SNRs are more than the 
results calculated from the two-fluid model. As for the en- 
ergy spectrum, the CR spectrum in this model shows a ba- 
sic power-law spectrum in the total energy range with a con- 
cave curve at some energy range because of the precursor 
structure. More detailed new st udies of the kinetic model can 
be referred to in these papers fKang. Jones & Gieselerl 12002; 
iKang & Jonesl 120071: lAmenomori et al.l 120081: iLi et al.l |2009; 
IZirakashvili & Aharonia n 2010). 

In an effort to follo w and extend the previo us dynamical 
Monte Carlo simulation (IKnerr. Jokipii & Ellisonll 1996k . we in- 
dependently developed a simulation code based on the Matlab 
platform using multiple scattering laws. Our multiple scatter- 
ing angular distributions consist of three Gaussian distributions 
and one isotropic distribution for the scattering angles during 
the scattering process. The aim of the isotropic scattering angu- 
lar distribution is to check the dynamical Monte Carlo method 
independently. Besides this, we want to know how the Gaussian 
distributions affect the scattering angular distribution function 
and the shock wave's evolution and propagation; even more, we 
expect to find the relationships between the multiple scattering 
law and the shock compression ratio. To validate the multiple 
scattering angular distributions, we followed the parallel-plane 
collisionless shock and the particle's acceleration using the same 
parameters and data as from Earth's bow shock, which was used 
in previous dynamical Monte Carlo simulation. 

In Section^ we outline the motivations for performing these 
four cases with four different scattering angular distributions. In 
Section [3J the specific simulation techniques are described. In 
Section |U we present the shock simulation results and different 
cases with four assumptions for scattering angle distributions. 
Section [5] includes a summary and the conclusions. 

2. The model 

The dynamical Monte Carlo si mulation has been developed by 
IKnerr. Jokipii & E llison (1996) to study Earth's bow shocks. It 



gives good results for the higher than lMeV cutoff in energy par- 
ticles and the power-law energy tail in the energy spectra. The 
dynamical Monte Carlo simulation method uses the prescribed 
scattering law instead of the complex electromagnetic field cal- 
culation like in the hybrid model. In addition, the dynamical 
Monte Carlo simulation need not assume the CR's injection rate 
and the associated diffusive coefficient as do the two-fluid and 
kinetic models. For the above reasons, we consider that develop- 
ing a simulation code by following the previous dynamical sim- 
ulation is necessary. Although the previous results successfully 
agree with observed data, the authors mention that their results 
show that the total compression ratio of the shock is more than 
4, which should be less than the ratio of the standard value for 
a nonrelativistic shock (Pelletier 2001). The Rankine-Hugoniot 
(RH) jump conditions allow deriving the relation of the compres- 
sion ratio with the Mach number: r = (y a + l)/(y a - 1 +2/M 2 ), for 
a nonrelativistic shock, the adiabatic index y a = 5/3 , if the Mach 
number M » 1 , then the maximum compression ratio should be 
4. To validate these consistent results from the previous model 
and extend this study to find what might be responsible for the 
shock compression ratio, we extend the previous isotropic scat- 
tering angular law by including an anisotropic scattering angu- 
lar law. This prescribed multiple scattering law consists of an 
isotropic scattering angular distribution and an anisotropic scat- 
tering angular distribution. The scattering angles consist of two 
variables of pitch angle and azimuthal angle. Once a particle has 
a collision with the massive scattering centers, its pitch angle 
becomes 6' '=6+66, and the azimuthal angle becomes cp'=cp+6(p, 
where 86 is the variation in the pitch angle 6, and Sep is the varia- 
tion in the azimuthal angle 0. The pitch angles 6 and 6 f are both 
in the range < 6, 6' < n, and azimuthal angles and 0' are both 
in the range < 0, 0' < 2n on the unit sphere. The variation in 
the pitch angle 66 and azimuthal angle Sep are composed of the 
scattering angle, and its anisotropic character is described by the 
Gaussian function f(66, Sep). 

Under the multiple scattering angular distribution law, four 
cases are calculated with three Gaussian distributions and one 
isotropic random distribution for the scattering angles. Here, the 
sign o~ is used to represent the standard deviation of the Gaussian 
function, and the sign ji is used to represent the statistical average 
or expected value of Gaussian function for the scattering angles 
(66, 6cp). We catalog the four cases as follows. 

(1) Case A: the scattering angles (66, 6cp) are distributed with 
a standard deviation <x = n/4 and an average value \i = 0. 

(2) Case B: the scattering angles (66, 6cp) are distributed with 
a standard deviation <x = n/2 and an average value fi = 0. 

(3) Case C: the scattering angles (66, 6(p) are distributed with 
a standard deviation <x = n and an average value fi = 0. 

(4) Case D: the scattering angles (66, 6cp) are distributed with 
a standard deviation cr = oo and an average value \i = 0, with 66 
varying from -n/2 to n/2, and 6cp varying from -n to n isotrop- 
ically. 

We performed four simulations according to the four differ- 
ent assumptions of the scattering angular distributions algorithm. 
We also assumed the scattering time (i.e., the mean time between 
two scattering events) is the same constant in the four cases as 
in the previous model. The idea that such a simple law can be 
used to descr ibe the entire scattering process was postulated by 
lEichled (1 19791). b ased on the two- stream instability work done 
by iParkerl (1 19611) . Put simply, it is assumed that the turbulence 
generated by both energetic particles streaming in front of the 
shock and by thermal particles produces nearly elastic scattering 
for particles for all energies in diffusive shocks. 
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3. Description of the method 

For simulating the total properties of shocks as they evolve from 
formation to a final steady state as energy increases via Fermi 
acceleration, we used the dynamical Monte Carlo model which 
employed the PIC techniques. Because there is no assumption 
of the injection rate or transparency function in PIC techniques, 
the shock-heated downstream ions can freely scatter back across 
the subshock into the upstream region without being thermal- 
ized, and the superthermal particles are produced in the thermal 
background self-consistently. In addition, unlike the hybrid sim- 
ulation, there is no complicated electromagnetic field calculation 
for individual particles, because it is replaced by the prescribed 
scattering law (Ellis on et.al.| [l993). To reproduce these acceler- 
ation and scattering processes, a similar simulation box and the 
same parameters (see Appendix |A| as the previous dynamical 
Monte Carlo method are used in these new codes. As described 
in the previous simulation (see Figure lAJl , the particles with an 
initial bulk velocity Uq and a Maxwellian thermal velocity Vl in 
their local frame are moving along a parallel magnetic field Bo m 
a one-dimensional box. To maintain a continuous flux- weighted 
flow upstream, a new particle fluid with the same density needs 
to be injected into the simulation box from the left boundary. For 
a shock initialization, the reflecting wall on the right boundary is 
used to reflect the incoming particles, and forms a piston shock. 
The model also includes the escape of energetic particles at the 
upstream free escape boundary (FEB). The FEB phenomenolog- 
ically models a finite shock size or the lack of sufficient scatter- 
ing far upstream to turn particles around. Once ions cross the 
FEB, they are assumed to decoup le from the shock sy stem, and 
are taken as the energy losses (I Jones & Ellison! fl99lf) . The size 
of the foreshock region (the distance from the shock front to the 
FEB) thus sets a limit on the maximum energy a particle can 
obtain. 

As shown in Figure IA.U one particle's box frame velocity 
V is a total velocity, which is composed of the local thermal 
velocity Vl and the bulk fluid speed U (i.e. V - Vl + U, for 
upstream U = Uq, for downstream U = 0). After one particle 
arrives in the downstream region, its kinetic energy is converted 
into random thermal energy by dissipation processes. With the 
development of these many processes, the bulk fluid speed of 
downstream flow becomes zero, and the length of downstream 
region is extended dynamically. 

As listed in Table I A. ll all of the specific parameters are used 
in our simulations, considering PIC techniques. The total length 
(X max = 300) is divided into the number of grids (n x = 600) 
with a grid length (Ax = 1/2). Initial grid density of the parti- 
cles (no = 650) is set. The total time (T max = 2400) is divided 
into the number of time steps (N t - 72000) with an increment 
of time (dt = 1/30). In summary, these new codes consist of 
the following three substeps like the previous simulation, except 
for the third substep employing the extended multiple scattering 
laws. 

(i) Individual particles move. Particles with their velocities 
move along the one-dimensional x axis: 



X[ = X{- 1 + (Vj k dt, te[l, t max ], k e [l,k max ], 
where 

(V x ) k = (V Lx ) k + (U k ) x , 



(U k ) x 



-f(v x ), 



(1) 



(2) 



(3) 



here r majc =72000, k max =600, and "k" represents the index of the 
computational grid, (Uk) x represents the bulk fluid speed of the 
computational grid along to the x direction, and the value of Uk 
is obtained from substep (ii). Since we are simulating a diffusive 
shock based on a one-dimensional parallel magnetic field, the 
fluid quantities only vary in the x direction. 

(ii) Mass collection. Summation of particle masses and ve- 
locities are calculated at the center of each computational grid: 



"A 

P k = Yj m p( v ^ k = 1.2,...** 



(=1 



Uk 



- V(n-,* = i,2,...*„ 



(4) 



(5) 



where is the number density of particles in the "k" grid, rep- 
resenting the mass of the computational grid. Here, Pk is the 
total momentum of the protons in the "k" grid, m p is the mass of 
an individual proton, and Uk is the average bulk fluid speed of 
the grid (i.e. the velocity of the scattering center). The collected 
grid-based mass and momentum densities will directly decide 
the velocity of the scattering center Uk. The particle's total ve- 
locity V in the box frame is decided by Equation 0. Once the 
value of Uk becomes zero, the shock front is decided by the po- 
sition of the corresponding grid, and it means that the shock is 
formed and the length of the downstream region is extended dy- 
namically. Similarly, if the value of Uk is between Uq and zero, 
it means that the foreshock region or precursor (i.e. between the 
FEB and the shock front) is formed by the "back pressure." The 
FEB and the shock front both dynamically move away from the 
reflective wall with a shock velocity v^. 

(iii) Applying multiple scattering laws. A certain fraction of 
the particles are chosen to scatter the background scattering cen- 
ter with their corresponding scattering angles according to the 
prescribed scattering angular distributions. The average number 
of scattering events occurring in an increment of time dt depends 
on the scattering time scale r, and the scattering rate is presented 
by 



R s = dt/r, 



(6) 



where R s is the probability of the scattering events occurring in 
an increment of time. The candidates with their local velocities 
and scattering angles scatter off the grid-based scattering centers. 
These individual particles do not change their routes until they 
are selected to scatter once again. So the particle's mean free 
path is proportional to the local thermal velocities in the local 
frame with 



AocV L , 

for simplicity, we take its formula as 
A = V L -r. 



(7) 



(8) 



For the individual protons, the grid-based scattering center can 
be seen as a sum of individual momenta. So these scattering pro- 
cesses can be taken as the elastic collisions. In an increment of 
time, once all of the candidates complete these elastic collisions, 
the momentum of the grid-based scattering center is changed. In 
turn, the momentum of the grid-based scattering center will af- 
fect the momenta of the individual particles in their correspond- 
ing grid in the next increment time. One complete time step con- 
sists of the above three substeps. The total simulation temporally 
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evolves forward by repeating this time step sequence. To calcu- 
late the scattering processes accurately and produce an exponen- 
tial mean free path distribution, the time step should be less than 
the scattering time (i.e. dt < r). 

The presented multiple scattering law simulations are devel- 
oped on the Matlab platform. Any one of the four cases can oc- 
cupy the CPU time for about seven weeks on a 3. 4GHZ (MF) 
CPU per core. To speed up the running programs, the parallel al- 
gorithm should be used on a high performance computer (HPC). 

4. Results & discussions 

We present all of the shock profiles for the shock simulations of 
the four cases in Figure [T] and we present all aspects of simula- 
tion results including the density and velocity profiles, compres- 
sion ratios, analysis of the heating and acceleration, and energy 
spectrum, as well as the correlations between the shock com- 
pressions with the energy spectral index. For the convenience of 
comparison and discussion, we list the specific calculated items 
in Table [T] Here, cr is the standard deviation of the scattering an- 
gular distribution function f(66, Sep), and fi is the average value 
of scattering angles (66, 6cp). The subshock compression ratios 
r su b are calculated from the velocity profiles in the same shock 
frame reference. The compression ratios r u and r p are calculated 
from velocity profiles and density profiles, respectively. The to- 
tal energy spectral index T tot and the subshock' s energy spectral 
index T su b are calculated from compression ratios r u and r su b, 
respectively. The last two rows are shown as scaled values. 

As shown in Figure [T] the present isotropic Case D largely 
appears similar to the results from prev ious dynamical simula- 
tions by lKnerr. Jokipii & Ellison! (Il996l) . In addition, all aspects 
of the shock wave structure, density and velocity profiles, com- 
pression ratio and energy spectra present in isotropic Case D 
also give similar results to the previous outcome. The specific 
results of the present isotropic Case D are shown in the fifth col- 
umn of Tableffl v^=-0.0733, X sh =124, FEB=34, r^=3.9444, 
r M =5.0909, r p ~ r u with a difference of 0.18%, r^=1.0094, 
r^=0.8667, | V^|=98.46km/s, and £ CMto// =4.01MeV. As a com- 
parison, the corresponding results of the previous simulation 
are also given: v^=-0.0720, X s h= 127.5, FEB=31.5, r su b=3.20, 
r M =5.20, r p ~ r u with a difference of 2.3% , r^=1.1818, 
r to ,=0.8571, \V sh \=96.9 km/s, and E cutoff ~4.00MeV. Among 
comparable results, a slightly larger differences in the values of 
the r su b and T su b between the two isotropic simulations would 
be distributed between different sizes of the subshock region, 
decided differently. As seen from the comparison of the results 
coming from two independent simulation codes, the present sim- 
ulation code successfully produced good agreement in the re- 
sults with those in previous dynamical Monte Carlo simulation. 
Therefore, the present simulation code is based on the Matlab 
platform without using a supercomputer that can independently 
validate the previous dynamical simulation method using a com- 
pletely different code for that supercomputer. Next, we offer a 
series of discussions about the different cases considering the 
specific aspects of the simulation results for diffusive shock. 

4.1. Shock Profiles 

Figure Q] shows time sequences of the density profiles of four 
cases. In each plot, a shock forms and moves away from the re- 
flective wall, and the dashed line represents the FEB position 
with the time in each case parallel to the shock front position. 
We can see that both the shock position and the FEB position 
are moving with a virtually constant velocity from the beginning 




Fig. 1. Four cases of density profiles for the entire simulation 
box vs the time. The dashed line represents the position of the 
FEB in each plot. 

of the simulation to the end of the simulation (i.e. T max =2400) 
in each case. Simultaneously, as far as the positions of the FEB 
are concerned, we can see that the FEB position at the end of the 
simulation is significantly different in four different cases. As for 
the average density fluctuation in the downstream region, there 
are also apparent changes in different cases, and case A has the 
slowest shock propagation speed among these four cases. Case 
D has the lowest average density profile in the downstream re- 
gion among these four cases. Because from Cases A to D the 
only difference is the prescribed scattering angular distribution, 
we conclude that these differences of the results for shock prop- 
agation speed and density profiles are decided by the standard 
deviation value cr of the scattering angular distribution. 

Figure |2] shows four cases of density and velocity profiles at 
the end of the simulations. From Cases A to D, the position of 
the FEB approaches zero as the value of the standard deviation cr 
increases. The effects of the accelerated particles are clearly seen 
in the upstream smoothing of the velocity profiles in each case. 
In the simulations, when high-energy particles cross the shock 
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Table 1. Results of calculation with an initial energy of £o=1.3105keV. 
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front and diffuse upstream, they contribute negatively to the ve- 
locity profile. This reduces the grid-based velocity in the zones 
upstream of the shock, which in turn affects particles that are 
scattering in that region. In fact, the accelerated particles slow 
and heat the incoming flow and smooth the shock transition by 
their "back reaction." As is obvious to see from the velocity and 
density plots, different scattering angular distributions produce 
different effects on the shock wave evolution. For the examples 
presented here, we consider that a difference of approximately 
40.93% of the shock velocity is contributed by the scattering an- 
gle distribution. 

4.2. Compression ratios 

Here, we compare the compression ratios calculated from the 
velocity profiles with those from the density relationships. First, 
the value of the total compression ratio can simply be calculated 
from the formula 

r u = ui/u 2 (9) 

where u\ = uo + \v s h\ , u 2 = \v s h\, and u\(u2) is the upstream 
(downstream) velocity in the shock frame. The shock velocity at 
the end of the simulation (T max ) can be derived from the formula 

Vsh — (X max — X s h)/T max (10) 

where X max =300, T max =2400, and X s h is the position of the 
shock at the end of the simulation (see Table []])• The specific 
calculated results are shown in Table [TJ 

But in terms of the specific shock structure as seen in Figure 
[3 an accurate subshock compression ratio calculation should be 
more complicated. In any one of the cases in Figure [3] (plotted in 
the box reference frame), we show the specific aspects of a shock 
modified by an energetic particle population whose mean-free- 
path is an increasing function of momentum. The shock structure 
in each plot consists of three main parts: precursor, subshock, 
and downstream. The smooth precursor is on the largest length 
scale between the FEB and near shock position X s h, where the 
fluid speed gradually decreases from value Uq to v su b. The size 
of this precursor is almost the mean-free-path length of the max- 
imum energy accelerated particles. One of the smallest scales is 
the subshock region with a sharp deflection of the fluid speed de- 
creasing from v su t to vtox = 0. The downstream region changes 



after the fluid speed becomes v\> ox - by microphysical dissi- 
pation processes. The gas subshock is just an ordinary discon- 
tinuous classical shock e mbedded in the comparably larger scale 
energetic particle shock (Berezhko & Ellison 1999). The value 
of v su b is determined by a sharp deflection of smooth curves in 
velocity profiles near the shock front, and the value of the sub- 
shock velocity increases from cases A, B, and C to Case D (i.e. 
(v su b)A < (v S ubh < (v S ub)c < iy sub ) D ). All of the velocity profiles 
are based on the box frame. That value of the box frame's ve- 
locity is zero (v^ox = 0) in all cases. The subshock compression 
ratio r su b is calculated from the formula r su b - (v su b + |v^|)/|v^| 
. For the sake of the comparison of the values of r su b in differ- 
ent cases, the subshock compression ratios are calculated in the 
same shock frame reference, and the calculated results of r su b are 
shown in Table [T] 

We then have the calculations of the compression ratio from 
the density relationships between the upstream and downstream 
flows: 

r p =p2/pi (11) 

where p\ = no is the upstream density, and p2 is the downstream 
density. This value is presented by which is the average value of 
rik over the downstream region. 

V^max ^shJ ^ 

where is the number density of particles in the "k" grid, k s h = 
x s hldx is the grid index of the shock at the end of the simulation 
(T max = 2400), and k max = 600 is the grid index of the X max . 

Figure IH shows the complete density plots of the four cases 
at the end of the simulation. The value of the upstream density p\ 
is the same constant value, which is equal to the initial density no 
in each case. The value of the downstream density p2 decreases 
from cases A, B, and C to Case D (i.e. (p2)a > (Pi)b > (p2)c > 
(P2)d)- Similarly, the detailed calculation results of the compres- 
sion ratios r p are listed in Table [T] As listed in Table [T] the values 
of the subshock compression ratios, r^=2.5421, r su b =3. 0207, 
r SM fc=3.3903, and r^=3.9444, corresponding to cases A, B, C, 
and D, respectively, are all lower than the standard value of 
r = 4. Unfortunately, the values of total compression ratio r u 
and r p in each case are both higher than the standard value of 
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Case A (a=rc/4,n=0) 




> 4000 
1 



D 50 100 150 200 250 3C 
Case B (a=7t/2,|i=0) 


ww^v^" 







250 300 





Fig. 2. The velocity and density profiles (in the box frame) vs 
the position at the end of the simulation (T max =2400) in the four 
cases. The vertical dashed line in each plot represents the posi- 
tion of the FEB, from case A to case D. 



r - 4. But Knerr. Jokipii & Ellison! (Il996l) consider that, if en- 
ergy is lost from the system (e.g. by particles escaping via FEB), 
it is possible to produce a shock with a total compression ratio 
that is higher than the standard value predicted by the Rankine- 
Hugoniot (RH) conditions. We have examined the mass and en- 
ergy losses via the FEB in each case. The results definitely show 
that the case with more energy losses would produce a higher 
total compression ratio than those in the case with less energy 
loss. Consequently, we consider that the energy loss rates would 
be affected by the prescribed scattering law. In any case, the en- 
ergy losses are always an important and interesting problem in 
the nonlinear diffusive shock acceleration theory, so we will per- 
form more precise research focusing on these problems in later 
papers. In addition, although the values of r u are correspondingly 
slightly higher than the value of r p in each case, all these differ- 
ences are less than 3%, and the specific difference in each case is 
2.9%, 1.5%, 1.3% and 0.18% corresponding to cases A, B, C and 
D. As seen from Figs. [3] and ffl the value of the total compression 
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Fig. 3. Velocity profiles in the shock region at the end of the simula- 
tion (r maJC =2400) in the four cases; the vertical solid and dotted lines 
indicate the shock front and FEB in each plot, respectively; the hori- 
zontal solid, dotted, dot-dashed and dashed lines show the values of the 
shock velocity v sh , velocity of box frame v box , subshock velocity v sub 
and initial bulk velocity U , respectively. Two vertical bars in each plot 
represent the two deflections of velocity, the upper bar represents the 
part of the shock precursor and the lower one represents the subshock. 
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Fig. 4. Density profiles in an entire simulation box at the end of 
the simulation (T max =2400) in the four cases. The vertical solid 
line is located in the position of the shock front, the upper hor- 
izontal dot-dashed line represents the value of the downstream 
density p2, and the lower horizontal dashed line indicates the 
value of the upstream density p\ in each plot. 



ratio r tot , determined from the velocity profiles, is more consis- 
tent with the density profiles in each case (i.e. the total compres- 
sion ratios r tot in all cases are satisfied by the Rankine-Hugoniot 
(RH) conditions : u\/u2 = pilp\). Therefore, it is not difficult 
for us to conclude that the total compression ratios r u and r p de- 
crease as the value of the standard deviation <x of the scattering 
angular distribution increases, but the subshock compression ra- 
tio r su t increases as the value of the standard deviation cr of the 
scattering angular distribution increases. 

4.3. Heating & acceleration 

Here we contrast between two important aspects of the heating 
and acceleration processes in the diffusive shock acceleration. 
Figure [5] shows the particles' scatter plots in four cases at the 
end of the simulations (T max =2400). In each case, particles with 
local velocity scatter into the simulation box's position. A large 
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Fig. 5. The scatter plots of the particle's thermal velocities in 
the local frame vs its position at the end of the simulations 
(Tmax=2400), and the vertical dashed line and solid line in each 
plot indicate the approximate position of the FEB and shock 
front, respectively. Only the ratio of 1/12 of the total number 
of particles are plotted. 



amount of particles that do not get injected into the Fermi accel- 
eration mechanism and that have lower thermal velocities stay 
in the downstream region, and a few particles with higher en- 
ergy via multiple shock encounters can move far away from the 
shock front and even escape from the FEB. From Cases A to D, 
more and more particles are injected into the Fermi acceleration 
mechanism, and they gain greater and greater maximum energy. 
Obviously, the maximum thermal velocity in Case D would be 
several times that of the ones in Case A. We supposed that this 
difference is mainly contributed by the scattering angular dis- 
tribution function f(S0, S(p). In short, the majority of particles, 
which flow toward the shock, cross the shock only once, and (af- 
ter scattering) remain fairly stationary in the downstream region, 
which would consist of the "heated" elements, and a few high- 
energy particles represent the "power-law" part of the simulated 
particles flows. Actually, the "back pressure" from the acceler- 
ated particles via their back reaction reduces the incoming fluid 
speed, leading to a smoothed precursor heating. Therefore, an 
anisotropic scattering angular distribution in the shocks dom- 
inate the "gas heating" process in the simulated plasma, and 
isotropic scattering angular distribution in the shocks play an 
important role in nonlinear acceleration of the energetic parti- 
cles via the Fermi mechanism, and they dominate the "precursor 
heating". 



Fig. 6. The final energy spectrum in the four cases and the initial 
energy spectrum plots. The thick solid line with a narrow peak at 
Eq =1.3105keV represents the initial Maxwell energy distribu- 
tion. The solid, dashed, dot-dashed, and dotted extended curves 
indicate the simulated particles' energy spectral distribution, av- 
eraged over the entire downstream region, at the end of the sim- 
ulations (T max =2400), corresponding to Cases A, B, C, and D, 
respectively. Most particles cross the shock only once, produc- 
ing the large broad peak centered at E^ ~ 0.05keV, E B ~0. IkeV, 
E c ~ 0.15keV, and E D ~0.20keV in Cases A, B, C, and D, re- 
spectively. However, some particles gain enough energy via the 
Fermi acceleration mechanism to produce the "power-law" tail 
in the energy spectrum with the cutoff at Ea= 1 . 10 MeV, £"5=2.41 
MeV, £ c =2.98 MeV, and E D =4.01 MeV corresponding to Cases 
A, B, C, and D. All these energy spectra are plotted in the same 
shock frame. 



4.4. Energy spectra 

Spectra calculated in the shock frame from the initial and fi- 
nal particle (proton) energy distributions in all cases are shown 
in Figure [6] The energy units in this plot are derived from the 
scaling parameters presented in Table IA.1I (see Appendix |A|). 
Initially, all particles move toward the wall with a certain thermal 
spread in energy. A narrow peak at E=1.3105keV represents the 
initial Maxwell energy distribution. The four extended curves 
indicate the simulated particle energy spectral distribution, aver- 
aged over the entire downstream region, at the end of the sim- 
ulations, corresponding to the four cases, respectively. The ma- 
jority of the particles cross the shock only once, producing an 
expanded energy spectrum with a central peak at E^ ~ 0.05keV, 
E B -O.lkeV, E c ~ 0.15keV, and E D ~0.20keV in Cases A, B, 
C, and D, respectively. However, as is shown in Figure [6l the 
minority of the particles gain enough energy via the Fermi ac- 
celeration mechanism to produce the "power-law" tail in the en- 
ergy spectrum with the cutoff at E A =1A0 MeV, E B -2A\ MeV, 
£ c =2.98 MeV and E D =4.01 MeV corresponding to Cases A, 
B, C and D, respectively. For more details about the calculated 
results, see Table [TJ It is evident from Figure [6] that the values 
of the central peak of the extended energy spectra in the four 
cases are far from the initial energy peak in their respective or- 
der. This means the values of the central peak in each case in- 
crease as the value of the standard deviation of the scattering 
angular distribution increases, and each extended curve shows 
a harder power-law slope in its high-energy tail as the expand 
energy range increases, respectively. Therefore, we can see that 
the case of applying an anisotropic scattering angular distribu- 
tion function will produce a slightly softer energy spectrum, and 
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the case of applying an isotropic scattering angular distribution 
will produce a slightly harder energy spectrum. 

4.5. Spectral index & compression ratios 

Usually, we could predict the power-law energy spectral index 
from diffusive shock acceleration theory: 



dJ/dE oc E~ Y 



(13) 



where dJ/dE is the energy flux and Y is the energy spectral in- 
dex, and 



r = (r + 2)/(2x(r-l)). 



(14) 



According to Equation[l4l we substituted the values of the com- 
pression ratio r with two group values of r tot and r su b obtained 
in each case. Then, the two group energy spectral indices T tot 
and T su b in each case are calculated. Two groups' spectral in- 
dex values are listed in Table Q] as Ta= 0.7167, T B =0.7677, 
T c =0.8083, and T D =0.8667 in the total group Y toU and T A = 
1.4727, T B =1.2423, T c =1.1275, and T D =1.0094 in subshock 
group T su b, corresponding to the cases A, B, C, and D, respec- 
tively. As shown in Figure[71 from Cases A to D, all of the values 
of the subshock' s energy spectral index T su b > 1 and show that 
a slightly harder power-law slope in the respective order. From 
Cases A to D, all of the values of the total energy spectral index 
Y tot < 1 show a slightly decreasingly deviation from the power- 
law slope in the respective energy spectrum. 



Spectral Index with Multiple Scattering Angular Distributions 




1/4 1/2 3/4 1 2 3 

a (The Deviation of Gaussian Distribution) 



Fig. 7. The correlation of the deviation value of the Gaussian 
distribution vs the energy spectral index. The triangles represent 
the total energy spectral index in each case. The circles indicate 
the subshock's energy spectral index in each case. From Cases 
A to D, all of the values of the subshock's energy spectral index 
^sub > 1 show a slightly harder power-law slope in the respective 
order. From Cases A to D, all of the values of the total energy 
spectral index T tot < 1 show a slightly decreasing deviation from 
the power-law slope in the respective energy spectrum. 



5. Summary and conclusions 

We followed the previous dynamical Monte Carlo simulation by 
using a new code based on the Matlab platform independently, 
and presented the same results as the outcome from the previous 
simulation. In addition, we successfully extended the simulation 
include the multiple scattering angular distributions using these 
new codes to study the diffusive shock acceleration mechanism 
further. 



In conclusion, the comparison of the calculated results come 
from different extensive cases, we find that the total energy spec- 
tral index increases as the standard deviation value of the scat- 
tering angular distribution increases, but the subshock's energy 
spectral index decreases as the standard deviation value of the 
scattering angular distribution increases. In these multiple scat- 
tering angular distribution simulation cases, the prescribed scat- 
tering law dominates the shock structure and plays an impor- 
tant role in balancing whether the particles have more heating 
or more acceleration. In other words, the cases with anisotropic 
scattering distribution give the overall velocity-deflection pre- 
cursor sizes, which are larger than the isotropic case, and give a 
relatively greater "heating" effect or less "acceleration" effect on 
background flows than does the isotropic case. 

As a result, the shock compression ratio and the energy spec- 
tral index are both modified naturally by the prescribed scat- 
tering law. Specifically, the cases of applying an anisotropic 
scattering distribution function will produce a slightly softer 
subshock's energy spectral index, and the case of applying an 
isotropic scattering angular distribution will produce a slightly 
harder subshock's energy spectral index. Simultaneously, from 
the isotropic case to the anisotropic case, the total energy spec- 
trum shows an increasing deviation from the "power-law" distri- 
bution. 

In addition, although we find no case producing the total 
compression ratio which should be less than the standard value 
4 according to the Rankine-Hugoniot (RH) jump conditions, the 
fact is clear that the prescribed scattering angular distribution 
function would have an effect on the total compression ratio. If 
there is a suitably prescribed scattering law that leads to much 
less energy loss, it is possible to constrain the total compression 
ratio to be less than 4. 
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The Schematic Diagram of the Simulation Box 



FEB Shock 

(Free Escape Boundary) 



200 250 300 X 

Reflective Wall 



Downstream 



Fig. A.l. Shock is produced by supersonic flow toward the sta- 
tionary reflective wall to the right. Continuous inflow of new 
particles occurs at the left boundary. Inflow velocity is t/o, a 
particle's total velocity is V, the local frame particle velocity 
is Vl, and thermal velocity < Vl >= vq. The two circles rep- 
resent one typical particle in the upstream region and one in 
the downstream region, respectively. The vertical solid line rep- 
resents the shock front, the vertical dashed line represents the 
FEB, the velocity of the shock is V shock, size of the foreshock 
region is Xfeb = 90, the upstream flow velocity U = Uq, and 
the downstream flow velocity is U = 0. The magnetic field Bq 
and inflow velocity Up are bo th normal to the shock front (see 
iKnerr. Jokipii & Ellison|[l996h . 



Appendix A: Simulation box & parameters 

With respect to the validity and consistency of verifying the previous dynamic 
Monte Carlo simulation method, the present simulation program uses the same 
simulation box and identical parameters as the previous dynamical Monte Carlo 
simulation ( Knerr, Jokipii & Elliso n 1996). The schematic diagram of the simu- 
lation box is shown in Figure |AH and all the simulation parameters are listed in 
TableCQl 
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Table A.l. The parameters of the simulated cases 



Inflow velocity 


Uq—\). j 


4UJKm/s 


Thermal speed 


Vq— U. UZ 


ZD. 7KII1/S 


Scattering time 


t=0. 833 


0. 13s 


Box size 


X max =300 


10tf e 


Total time 


T max =2400 


6. 3minutes 


Time step size 


dt=l/30 


0. 0053s 


Number of zones 


nx=600 




Initial particles per cell 


^0=650 




FEB distance 


X feb =90 


3R e 



Notes. Scaling used a box size = \0R e (where R e represents the Earth's radius) and the box frame inflow velocity uq- 403km/s. This implies 
the following scale factors for distance, velocity, and time: X sca / e = 10/^/300, v sca / e =403km s _1 /0.3, and t sca i e =x sca i e /v sca ie- Here, the Mach number 
M=11.6. Dimensionless or no rmalized numbers are used in the text to describe our simulations, except for specifically highlighted examples (see 
iKnerr. Jokipii & Ellisonll 19961) . 
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